Tree Density from LiDAR

library(bowen.biodiversity.webapp)
library(here)
#> here() starts at /home/jay/Programming_Projects/bowen.biodiversity.webapp
library(sf)
#> Linking to GEOS 3.11.1, GDAL 3.6.4, PROJ 9.1.1; sf_use_s2() is TRUE
library(leaflet)
library(leaflet.extras)
library(future)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

# install.packages('lasR', repos = 'https://r-lidar.r-universe.dev')
library(lasR)
#> 
#> Attaching package: 'lasR'
#> The following object is masked from 'package:dplyr':
#> 
#>     summarise
#> The following object is masked from 'package:future':
#> 
#>     sequential

Downloaded from LidarBC: https://lidar.gov.bc.ca/

Convert from .laz to .las

set_parallel_strategy(concurrent_files(4))
# Read las 
ctg <- readLAScatalog(here("temp/las"))

del = triangulate(filter = keep_ground())
norm = transform_with(del, "-")
write = write_las(here("temp/normalized/*_normalized.las"))
chm = rasterize(2, "max", ofile = here("temp/chm/chm.tif"))
lm = local_maximum_raster(chm, 5)
pipeline2 = reader() + del + norm + write + chm + lm

ans = exec(pipeline2, on = ctg)
ttops <- ans$local_maximum %>% st_transform(4326)
st_write(ttops, here("data-raw/ttops.gpkg"))

ttops <- st_read(here("data-raw/ttops.gpkg")) 
#> Reading layer `ttops' from data source 
#>   `/home/jay/Programming_Projects/bowen.biodiversity.webapp/data-raw/ttops.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 660834 features and 0 fields
#> Geometry type: POINT
#> Dimension:     XYZ
#> Bounding box:  xmin: -123.4483 ymin: 49.32728 xmax: -123.3 ymax: 49.42478
#> z_range:       zmin: 2 zmax: 1383.31
#> Geodetic CRS:  WGS 84
ttops_coord <- ttops %>% st_coordinates() %>% as.data.frame()
leaflet(ttops_coord) %>%      
  addProviderTiles(providers$CartoDB.Positron) %>%
  addHeatmap(lng=~X, lat=~Y, blur= 15, max = 65, radius = 30)